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INTRODUCTION 


This report summarizes the progress made during the period January 1 - 
June 30, 1993 under the NASA Grant NAG-3-768 titled "Numerical studiy of the 
effects of icing on viscous flow over wings". The work was carried out by Dr. L. N. 
Sankar, the principal investigator, with the assistance of two graduate research 
assistants, Mr. Ashok Bangalore and Mr. Napporn Phaengsook. Another student, 
Mr. Olympio Mello, not supported under this project, also contributed to the work 
reported. 

It may be recalled that the tasks to be performed during the current grant 
year are as follows: 

1 . Development of 3-D boundary layer methods for accurate estimates of the 
velocity field and surface heat transfer rates in the vicinity of the leading edge ice 
shape. 

2. Studies of the Effects of icing on 3-D High Lift System performance 

3. Continued improvement and validation of the 3-D Navier-Stokes solver 

Progress was made in these areas during the reporting period. In the 
appendix, abstracts of two papers submitted to the AIAA 1994 Aerospace 
Sciences Meeting may be found. These abstracts give details of the numerical 
formulation and several code validation studies. Here, only an overview of the 
work done to date is presented. 

Development of a 3-D Boundary Laver Methods: This work started out as a 
stand-alone 3-D steady, incompressible, boundary layer solver that will accept 
the surface pressure distribution from an external solver such as a panel method 
or a 3-D Navier-Stokes method, and will yield the velocity field and temperature 
field near the leading edge icing shape. 


It soon became apparent to us that a stand-alone 3-D boundary layer 
analysis, while computationally efficient, had a number of disadvantages. 
Specifically, 

a) 3-D boundary layer analyses break down when significant streamwise 
separation occurs. In the case of iced wing calculations separation can occur 
within the first 1 0% of the wing chord. 

b) steady boundary layer analyses can not model unsteady effects such as 
shear layer instabilities, vortex shedding etc. that have been reported by other 
researchers. 

c) The boundary layer analysis will have to be eventually coupled to an 
inviscid analysis, and an iterative method must be devised for full viscous-inviscid 
interaction. Such a fully interactive analysis can be very costly, and will break 
down if the wing stalls due to leading edge separation. 

Of course, a 3-D Navier-Stokes analysis is very costly and is also not an 
acceptable way of routinely solving 3-D viscous flows. 

For these reasons, an entirely new approach, based on some earlier work 
done by the present investigator was developed. In this approach, a two zone 
formulation is used. In the inner zone, the 3-D compressible Navier-Stokes 
equations are solved in their time dependent form. These equations do not have 
any singularity near the separation line, and can handle steady as well as 
unsteady flow phenomena. In the outer region, the 3-D unsteady compressible 
potential flow equation is solved. The inner and outer regions are tightly coupled 
to each other, so that a fully interactive analysis is achieved. There is no need to 
use a separate solver to specify the surface pressures, or to use viscous inviscid 
interaction iterations. Because the potential flow equation is solved in a fraction of 
the cost of the Navier-Stokes simulation, and because the Navier-Stokes zone is 
confined to a small region surrounding the wing, this hybrid formulation is 
extremely efficient, yielding results in 50% of the CPU time needed by our earlier 
Navier-Stokes simulations. 


The abstract of a paper by Mello et al. found in the Appendix describes the 
hybrid procedure in detail and gives some sample results for the new formulation. 

E ffects of icing on 3-D High Lift System Performance: To study this problem, we 
are making use of a multi-block version of the Georgia Tech 3-D solver. The 
multi-block solver is being set up so that an arbitrary number of flaps, slats and 
airfoils may be used in three dimensions. The flow field is divided into a number 
of zones or blocks in which the flow field is advanced by one time step. The block 
boundary may be a solid boundary, a block interface, a far-field boundary 
(inflow/outflow) or a symmetry plane. We have developed a fairly general way of 
setting up the blocks and specifying the block interface boundary conditions. 

At this writing, this solver is being debugged. We anticipate preliminary 
results from this solver to be available by September 1993. 

Continued Development and Validation of the 3-D Navier-Stokes Solver: It may 
be recalled that the Georgia Tech 3-D solver does a reasonable job of predicting 
surface pressure distributions over iced wings. The velocity field correlations are, 
however, not satisfactory. In an effort to improve the velocity field predictions of 
the 3-D Navier-Stokes solver, during the past nine months several improvements 
have been made to our solver. An approximate Riemann solver replaces the 
central difference scheme used in the earlier versions of our code, and a k-e 
turbulence model is available as an option in the solver, in addition to the 
Baldwin-Lomax base-line model. The appendix includes the abstract of a paper 
describing these enhancements, and application of this solver to iced wing 
configurations tested by Prof. Bragg. 
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SUMMARY 


A 3-D compressible Navier-Stokes solver has been developed and applied 
to 3-D viscous flow over clean and iced wings. This method uses a third order 
accurate finite volume scheme with flux difference splitting to model the inviscid 
fluxes, and second order accurate symmetric differences to model the viscous 
terms. The effects of turbulence are modeled using a k-e model. In the vicinity of 
the solid walls, the k and e values are modeled using Gorski's algebraic model. 

Sample results are presented for surface pressure distributions, for an 
untapered swept wing, made of NACA 0012 airfoil sections. The leading edge of 
these sections are modified using a simulated ice shape. Comparisons with the 
experimental data obtained by Prof. Bragg at the University of Illinois are given. 

The final paper will give these results, and additional fine grid calculations 
now in progress. 


INTRODUCTION 

The aerodynamic characteristics of lifting surfaces such as wings and 
horizontal tails may be dramatically altered by the accumulation of ice over the 
leading edge. Even a relatively small amount of ice accumulation may lead to 
dramatic rise in pressure drag, loss in lift, and in some cases, premature stall. 


Under the support of the NASA Lewis Research Center, a research effort 
has been underway at Georgia Tech on the effects of icing on the aerodynamic 
performance of lifting surfaces. Both 2-D and 3-D analyses have been 
developed, and calibrated [Ref. 1-3]. These analyses used relatively simple 
algorithms to solve the compressible Navier-Stokes equations. For example, 
standard second order accurate central differences are used to model the 
derivatives in the governing equations. A fourth order low-pass filter is used at 
every time step for removing any high frequency spatial oscillations. These 
analyses also used simple algebraic turbulence models such as the Baldwin- 
Lomax model to account for the effects of turbulence. Finally, for the sake of 
computational stability and efficiency, these methods used an implicit time- 
marching algorithm that permit the use of very large time steps. Reliable 
convergence to steady state in less than 1 CPU hour of computer time on a Cray 
Y/MP class of computer, on a 100,000 point grid. Reasonable correlation was 
observed between the calculations and the measured data obtained by Bragg 
and his coworkers [Ref. 4], 

As these analyses were applied to a number of iced wing calculations, two 
drawbacks of these approaches became apparent to the present investigators. 
First of all, the turbulence model used in these analyses is a simple algebraic 
eddy viscosity model, which uses the boundary layer thickness 6* (or an 
equivalent quantity such as the distance from the wall, y) as the length scale of 
turbulent eddies. The algebraic model used the boundary layer edge velocity U(x) 
(or an equivalent quantity y Icolmax) as a measure of the velocity scale of 
turbulent eddies. Such a simple representation of turbulence works well only for 
simple attached or mildly separated flows. In the case of iced wings, the flow field 
contains boundary layer over the solid surfaces, reversed flow regions, and a free 
shear layer coming off sharp corners and turns associated with the ice shape. 
The turbulent length and velocity scales for free shear layers are clearly different 
from the wall bounded flows. Use of a simple eddy viscosity model to such 
complex flows lead to inaccurate results in regions of separation. For example, 
the length of the separation bubble was often inaccurately predicted. Potapczuk 
has discussed the limitations of the Baldwin-Lomax model, and has 
recommended phenomenological fixes to this model [Ref. 5]. 
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The second drawback of the previous analyses is the irrecoverable 
filtering out and loss of information by the low-pass filters. These filters tend to 
smear out shock discontinuities, contact discontinuities, and shear layers. The 
rate at which the filtering is applied is proportional to lul+a , where lul is the 
magnitude of local flow, and a is the speed of sound. 

The present work improves upon the existing numerical analysis, 
presented in Ref. 3, and attempts to rectify the above two drawbacks. It uses a 
flux difference scheme, which dynamically and adaptively reduces the smearing 
of contact discontinuities, and vorticity waves, and leads to crisp resolution of free 
shear layers that arise in the flow. This method also uses a sophisticated k-e 
model to account for the effects of turbulence. In the vicinity of solid walls, this 
model uses a set of phenomenologically derived boundary conditions for k and e. 
The improved analysis allows transition to be detected automatically by the 
solver as regions of very low turbulent kinetic energy. 

The improved solver has been applied to clean and iced, swept wing 
configurations tested by Prof. Bragg and his coworkers at the University of Illinois 
[Ref. 4], 

The rest of this abstract is organized as follows. First, the Roe scheme is 
briefly outlined. Next, the k-e formulation is described. Several sample 
applications of the improved solver to iced wing analysis are presented next. This 
abstract concludes with a list of additional results to be found in the full paper. 

MATHEMATICAL AND NUMERICAL FORMULATION 

The 3-D compressible Navier-Stokes equations are given in a Cartesian 
coordinate system as 


dq dF dG dH dR as dT 

dt dx dy dz dx dy dz 

( 1 ) 

Here q is a vector of conserved flow properties; F, G and H are the inviscid 
flux vectors; R, S and T are viscous terms, which include the effects of turbulence 
through the eddy viscosity concept. 
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The calculations are carried out in a curvilinear body-fitted coordinate 
system (g l Ti£). |n this system, the governing equations may be written as 

dq dP dG dH_ dR dS df 

dt + + dr] + <?£ " d% dr] <?£ 

(2) 

where q etc. are related to their Cartesian counterparts by the metrics of 
transformation. 


The Roe Scheme: 


The objective of the time-marching algorithm is to advance the flow field at 
cells (or nodes (i,j,k) from a time level ’n' to the next time level 'n+1'. For this 
purpose, In the present work, equation (2) was discretized as follows: 


*n+l 

k ~~ Qi t j,k 


+ d g F;;i + 6 n G,% + 6 C HZ\ - + hK* 


(3) 


as 


where the operators etc. represent standard central differences such 


5,F 


F „„,- F ,-i.i 

A| 


( 4 ) 


The numerical flux vector T differs from the physical flux vector F in a 
manner specified by Roe [Ref. 6]. Let qi_ and qR be two estimates to the flow 
field vector q , to the left side and right side of the cell face (i+1/2,j,k). Then, 
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Kv2 - Hq ‘- )+ 2 F(q ^ ) + ^TW-'(q L - q R ) 
where. 


Qr “ 9i+l 
and 


(?i+l %. Qi+l ffi+l 


q L -q, + 3 M - + ~ — 

HL h ' 3 6 

(5) 

Here T is a matrix that contains the left eigenvectors of the matrix dF I dq, 
and A contains the eigenvalues of this matrix. These quantities must be 
computed using special "Roe" averages of qi_and qR. 

It may be shown that the above estimates of qi_ and qR lead to a spatially 
third order accurate formulation, on uniform grids and on mildly stretched grids. 


Equation (3) represents a system of non-linear algebraic equations for the 
flow properties vector q at time level 'n+1'. Furthermore, each node (i,j+1,k) is 
coupled to its nine neighbors. Thus, this system of equations are highly coupled. 
The classical Beam-Warming approximate factorization scheme was used solve 
these coupled nonlinear equations, in a manner similar to that described in 
Reference 1 . 

K-e Turbulence Model 


The starting point is the Boussnesq assumption for the turbulent stresses. 


-pUjUj - |i,( 


du, dllj 
dXj dXj 



The turbulent viscosity ^ has to modeled. This is achieved by modeling 

two physically conceivable quantities, the turbulent kinetic energy k, and the 
dissipation rate e. The basic high Reynolds number K-e model integrates the 
following transport equations for k and e. 


The turbulent production P, is given by: 


— du 


and the turbulent viscosity nj is related to k and e by : 

k 2 

Ht “ C ^ 

E 


The basic k-e model constants are as follows. 

C K -.09 
C t - 1.43 
C 2 - 1.92 
a k -1.0 
a E - 1.3 


The k-e Equations in Generalized Coordinates 


Before the equations are discretized, they have to be transformed to the 
generalized |, rj, C, coordinate system. In this coordinate system, the k-e equation 

set may be written as 
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Here U,V,W art the contravariant components of velocity and J is the 
Jacobian. 

In expanding the equations one reasonable assumption is made. In the 
diffusion terms, only the gradient of the quantities k and e in the £ direction are 
kept, the variation in the | and directions are neglected. This is justifiable as 
the maximum changes in the turbulent quantities occur only across the shear 
layer. The production term P is given by, 


C„k 2 P 

Je 


p f du} 3u + { du dv^ du | ( du ^ dw} 3u ^ 
l dx) dx \3y dx) dy \dz dx) dz 

I dv 3u^ 3v / 3v^ 3v ( dv 3w^ 3v 

\3x dy) dx \ dy) dy \3z 3y/ 3z 

( dw du} 3w ( dw dv} dw ^ f dw} dw 

1 3x dz) dx dy dz) dz V dz) dz 


The various derivatives with respect to x, y and z may be computed in terms of 
the generalized coordinates as follows, 


A A A r A 

dx " * x 31 +Tlx 3r| + ^ x 3t 
d d d y d 

3^“ ly 31 +Tly d^ +Cy ^ 

A t A A r A 

3z " §z 3l + 11z 3ri + ^ z 3t 
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The goal of the k-e solver is to compute k and e at a time level 'n+1 given 
these values at time level 'n'. In this formulation the "right hand side residual" is 
first calculated for the whole region . 

RHS = Diffusion - Convection + Production - Dissipation 

The convection terms were computed using two-point upwind schemes, as 
follows: 
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The diffusion terms were computed using central differences as follows: 

(t>')j.i/ 2 -(-V s ) jil , 2 '( k 4 


d f n t VrrVC k 1 _ J_ 
ar|[a e J ’•J Ati 


V °k 


'j+1'2 


d r ut Vt-VE 1_ 1 / 
dt, [o k J tj A^ 


V °k 




After the right hand side is computed, an approximate factorization 
scheme is used here to solve for the changes in k and e from tone time step. This 
requires the solution of the following matrix product equation. 
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[ 1 + t/A/<5 ? J 1 + VM5 n Jl + WAt6> 


- if r +1 -r 


e n+1 - e r 


- /A^KT/S)" 


Each of the three matrices (or operators) are invertible. 


To start the calculations, the k and e fields were initialized to small non- 
zero positive values. Near solid walls, the classical high Reynolds number k-e 
model breaks down. In the present work, near the solid wall, the k and e were 
assumed to have the following variation: 


k = Ay 3 


e = constant 


RESULTS AND DISCUSSION 

Navier-Stokes calculations were carried out using the flow solver 
developed in the previous section, for a swept untapered wing tested by Bragg et 
al. The airfoil sections are made of NACA 0012 airfoil sections. The wing tested 
had a chord length of 15 inches and a semispan of 37.25 inches. To study the 
effects of icing, the leading edge was modified with a simulated ice shape. This 
ice shape corresponds to that measured at the NASA Lewis Icing Research 
Tunnel, for a NACA 0012 airfoil, at a freestream velocity of 130 mph, angle of 
attack of 4 degrees, icing time of 4 minutes, volume median droplet diameter of 
20 microns, liquid water content of 2.1 grams per cubic meter and a temperature 
of 18 degrees. Bragg's measurements were done at zero sweep and at a 30 
degree sweep. 

The calculations used a 121 x 24 x 45 grid, with 91 points on the airfoil at 
each span station, 14 spanwise stations over the wing, and 45 points in the 
direction normal to the wing surface. Roughly a quarter of these points are in the 
boundary layer region. Thus, the boundary layer is only moderately resolved in 
these preliminary calculations. 
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Figure 1 shows the surface pressure distribution for the clean wing at 8 
degrees angle of attack at 5 spanwise locations: 27%, 42%, 56%, 72% and 89%. 
In general, the comparisons are in good agreement with the measured data. The 
suction peak near the leading edge, and the rapid compression of the flow over 
the upper surface were well predicted. No angle of attack corrections to account 
for the effects of the wind tunnel blockage were made. The flow solver had no 
user adjustable constants , and was run using the standard k-e constants 

described in the previous sections. 

Figure 2 shows the calculations and measurements for the iced wing. In 
this case, the flow separates over the leading edge ice shape, and forms a 
plateau. The Roe k-e is seen to perform satisfactorily at the inboard stations. 
Near the wing tip, the agreement is poor. At this writing, we are repeating this 
particular analysis with a clustered grid near the tip region to see if an improved 
resolution of the tip vortex will lead to an improved prediction of the surface 
pressures at the outboard station. 

Figure 3 shows some preliminary results for the velocity field at selected 
spanwise and chordwise locations, for the iced swept wing. Both the k-e results 
and the original algebraic model results are shown. The computed qualities 
plotted are the u- component (chordwise) of velocity, while the measured data 
are tangential to to the airfoil/wing surface. As a result, a discrepancy exists at 
the boundary layer edge. 

ADDITIONAL RESULTS TO BE PRESENTED IN THE FULL PAPER 

a) The full paper will contain results for the swept wing considered above 
and the unswept wing, on a denser grid. We are increasing the number of nodes 
in the boundary layer region, and in the spanwise direction, particularly near the 
tip. The final paper will contain these fine grid calculations. 

b) A wealth of LDV data is available in the Ph. D. dissertation of 
Khodadoust [Ref. 7]. The computed velocity profile at selected stations over the 
clean and iced wing will be compared with the measured data in the full paper. 
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CONCLUDING REMARKS 


An existing compressible Navier-Stokes solver for clean and iced wing 
analysis has been upgraded through the use of state of the art upwind difference 
schemes, and a two equation k-e model. Encouraging preliminary results have 
been obtained for swept, clean and iced wing configurations. The additional 
overhead associated with the Roe scheme and the k-e model has been partially 
offset by improvements to the algorithm elsewhere. A careful systematic 
validation of the improved flow solver is now underway, prior to its adaptation for 
iced wing performance analysis. 
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SUMMARY 

A hybrid method for computing compressible viscous flows is 
presented. This method divides the computational domain into two 
zones. In the outer zone, the unsteady full-potential equation (FPE) is 
solved. In the inner zone, the Navier-Stokes equations are solved. 
The two zones are tightly coupled so that steady and unsteady flows 
may be efficiently solved. The resulting CPU times are less than 50% 
of the required for a full-blown Navier-Stokes analysis. Sample 
applications of the method to an unswept iced wing at 4° angle of 
attack are presented. Surface pressures are in good agreement with 
the measurements obtained by Bragg et al. at the University of 
Illinois. The full paper will contain additional results for clean and 
iced swept wing configurations and will include surface pressure as 
well as velocity field comparisons. 


INTRODUCTION 

The increased need for aircraft to operate under adverse 
weather conditions has led to a growing effort in the performance 
evaluation of aircraft under icing conditions. Although the qualitative 
effects of icing have been well known to include lift decrease and 
drag increase with premature separation at the leading edge, there 
remains a need for reliable quantitative results in order to allow 
operation under icing conditions within an adequate safety margin. 
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The large number of parameters involved in describing the shape of 
ice accretion and its evolution in time make it impractical and 
expensive to build up the needed quantitative results solely from 
wind-tunnel studies, hence accurate numerical methods are clearly 
needed for this purpose. 

Current work in this area include the experimental efforts by 
Bragg et al. 1-3 , and the Navier-Stokes computations by Potapczuk 4 , 
Cebeci 5 and Sankar et al. 6 '9. In an ongoing research effort at Georgia 
Tech, a 3-D Navier-Stokes code has been successfully applied to iced 
wing configurations. While the Navier-Stokes analysis gives 
satisfactory results, it is CPU intensive. Since a typical iced wing 
performance would involve several parameters (a, Mo, Re, planform, 
ice accretion shape), the full Navier-Stokes code will be time- 
consuming. On the other hand, panel and FPE methods cannot 
accurately predict the loss in performance due to icing, even when 
coupled with a boundary-layer analysis, due to the occurrence of 
premature leading-edge separation and reversed flow. 

In order to improve the efficiency of Navier-Stokes 
computations, the present work was initiated. In this approach, a 
zonal Navier-Stokes/Full-Potential solver developed by Sankar et 
al. 10 and extended to rotors by Tsung and Sankar 11 is applied to an 
iced wing configuration. 

The Full-Potential solver used in the outer zone solves the 
unsteady compressible FPE in strong conservation form using the 
artificial compressibility concept and employing a strongly implicit 
procedure 12 . The Navier-Stokes solver used in the inner zone is 
identical to the one used in the previous iced wing efforts at Georgia 
Tech6-9. 


Historically, coupling potential flow to viscous flow via 
boundary-layer analysis has proved troublesome at the separation 
point and the recirculation region. Since we are computing the full 
Navier-Stokes equations in time-dependent form in the inner region, 
the above difficulties associated with boundary-layer methods are 
avoided. 

The re main ing of this abstract is organized as follows: First, the 
mathematical formulation of the Full-Potential and Navier-Stokes 
solvers are briefly given. Next, some sample applications of this 
method to the iced wing configuration tested by Bragg et al. 1 - 2 are 



given. The abstract concludes with a list of additional results to be 
presented in the full paper. 


MATHEMATICAL FORMULATION 


Full-Potenti al Formulation 

The 3-D unsteady compressible potential flow equation, in a 
body-fitted coordinate system, may be written in a strong 
conservation form as: 



( 1 ) 


where p is density and U, V and W are the contravariant components 
of velocity. 

Using the isentropic gas law and the energy equation, equation 
(1) may be written 12 as a second order hyperbolic partial differential 

equation for <t>: 



+ v 0„ x +W 



(2) 


At each point on the grid, the spatial derivatives are written 
using standard central difference formulas, which result in formal 
second-order accuracy in space. 

Because equation (2) is hyperbolic, it may be advanced in time 
using a stable marching scheme. Here, a strongly implicit procedure 

is used. To circumvent the nonlinearities, the coefficients p , a 2 , J , U , 
V and W appearing on the left side, and the density p appearing on 
the right side of equation (2) are computed at the time level 'n\ The 
remaining quantities in (2) are kept at the new time level 'n+1'. 

The time derivatives and the mixed time and space derivatives 
appearing on the left side of (2) are approximated using first order 
accurate finite difference formulas such as: 






H-1 


Ar 2 


(3a) 



(3b) 


A/A£ 

When the above discretizations are employed, at each time step 
a system of linear equations result for the quantity A <p = <p* + ' - <p\ This 
system of equations may formally be written in the form: 

[M] { A<J>} = {R} (4) 

The matrix M is a sparse matrix with twenty seven diagonals, 
reflecting the fact that a node (i,j,k) is linked to its 26 neighbors. On 
nearly orthogonal grids such as the one used in the present work, 
only diagonals coupling (i,j,k) to its six neighbors (i+l,j,k) , (i-l,j,k), 
(i,j-l,k), (i,j+l,k), (i,j,k-l) and (i,j,k+l) are significant, and one may 
neglect the other cross-coupling terms, which results in a seven- 
diagonal matrix. This seven-diagonal system is then recast into the 
LU form using a strongly implicit procedure. 

Navier-Stokes Formulation 

The vector form of the full Reynolds-averaged, 3-D Navier- 
Stokes equations based on an arbitrary curvilinear coordinate system 
can be written as: 

Q, + ( E - E v ) k + ( F - F v \ + ( G - Gy ) c = 0 (5) 

where Q. is the vector of unknown flow properties; E, F, G are the 
inviscid flux vectors; and E v , F v , Gv are the viscous flux vectors. 

The time derivative, Q*, of equation (5) is approximated using 
two-point backward difference at the new time level: 

Q.- tO"* 1 - « > (6, 

^ At 

where ‘n’ refers to the time level at which all quantities are known, 
and ‘n+1’ is the new time level. All spatial derivatives are 
approximated by standard second-order central differences and are 
represented by the differencing operators 5, e.g. E^= 5^ E. 

The spanwise derivative, F^, is evaluated explicitly at the old 
time level 'n\ but uses the ‘n+1’ values as soon as they become 



available. This semi-explicit treatment of the spanwise derivative 
enables the scheme to solve implicitly for all points at one spanwise 
station at a time. To eliminate any dependency the solution may 
have on the marching direction, the solver reverses the direction oi 
marching with every spanwise sweep. 


The viscous terms E v , F v and Gy are evaluated explicitly, using 
half-point central differencing, so that the computational stencil for 
the stress terms uses only three nodes in each of the three directions. 
Explicit treatment of the stress terms still permits the use of large 
time steps since the Reynolds numbers of interest here are fair y 

large. 

The time and space discretizations described above lead to a 
system of non-linear, block penta-diagonal matrix equations for the 
unknown Q. n+1 . The equation is then linearized using the Jacobian 
matrices A = 9E/aQ.and B = 3G/dQ.and approximately factored into a 
product of two block tri-diagonal matrix equations. 


[ I + At 6^ A ] [ I + At 5^ B ] AQ,= R 

= - At [ ( E - E v ) + 5 ^ ( F - F v ) + 6 ^ ( G - G v ) ] 


(7) 


The use of standard central differences to approximate the 
spatial derivatives can give rise to the growth of high frequency 
errors in the numerical solution with time. To control this growth, a 
set of 2nd/4th order non-linear, spectral radius based, explicit 
artificial dissipation terms are added to the discretized equations. A 
second order implicit dissipation is used to help the overall 
numerical stability of the scheme. 

A slightly modified version of the Baldwin-Lomax (B-L) 
algebraic turbulence model is used, where the maximum normal 
shear value is used instead of the wall shear stress because in the 
vicinity of separation points, the shear stress values approach zero at 

the wadi. 

Navier-Stok es/Full Potential Coupling 

A typical partitioning of the domain into an inner zone and an 
outer zone is illustrated in Fig. 1. The plane k=kMATCH corresponds to 
the interface between the inner zone and the outer zone. The Navier- 


Stokes solver is applied to all planes up to k=kMATCH- The FPE solver 
is applied between the planes k=kMATCH -1 and the outer boundary 
k-plane. Therefore, the two zones actually overlap, which allows 
specification of boundary conditions at the interface without 
extrapolation. 


RESULTS AND DISCUSSION 

The hybrid Navier-Stokes/Full-Potential Method has been 
applied to an iced unswept wing configuration, which has been 
experimentally studied by Bragg et al.1'3. The surface pressures were 
measured at five spanwise stations: 17%, 34%, 50%, 66% and 85%. 

The computational grid used in the present study was an 
algebraic C-grid with 141 x 19 x 41 grid points, with 121 points over 
the airfoil surface at each spanwise station. 14 spanwise stations 
were used along the wing, with 5 stations extending beyond the tip. 
The Navier-Stokes and Full-Potential solvers were interfaced at 
kMATCH=21. The Mach number was 0.12, the Reynolds number 1.7321 
million, and the angle of attack was 4°. The CPU time for this 
configuration was about 24 sec. per iteration on a Hewlett-Packard 
Apollo 700 workstation. The CPU time for the same configuration on 
Georgia Tech's Cray Y-MP/E was about 9.6 sec. per iteration. These 
times are about 50% of the times required for a full Navier-Stokes 
analysis. Computations were performed up to 1000 iterations, when 
the maximum residual and corrections were of order 10' 5 . 

The computed surface pressure distributions were linearly 
interpolated to allow comparisons at the spanwise stations where 
experimental data are available 1*2. The resulting correlations are 
shown in Fig. 2. A good agreement may be observed between the 
numerical and experimental results, especially inboard, where the 
suction peaks behind the ice accretion can be clearly observed, 
although the sharp variations are somewhat displaced with respect 
to the experimental data. Less agreement is found towards the tip 
region, probably due to the three-dimensional character of the flow 
in this region, which would require a finer grid resolution near the 
tip. There are also spurious peaks near the trailing edge at all the 
spanwise stations, indicating that the grid resolution is not adequate 
there. 


Overall, it is observed that our current results are basically 
identical to those obtained by a full Navier-Stokes code 6 * 7 * 9 , while 
consuming only about 50% of the CPU time. 


ADDITIONAL RESULTS TO BE FOUND IN THE FULL PAPER 

The final paper will include additional results for clean and 
iced wing configuration, including surface pressure data and 
correlation of velocity profiles with LDV flowfield measurements 
obtained by Bragg et al. 3 for both unswept and swept wing 
configurations. We will also address in more detail further 
developments of this hybrid Navier-Stokes/Full-Potential method. 
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Figures: 


Fig. 1: Partitioning of Computational Domain into Inner and Outer 
Zones. 




Fig. 2: Surface Pressure Distributions for Iced Unswept Wing at 4 
Angle of Attack. 



